// -*- Mode: C++; tab-width: 2; -*-
// vi: set ts=2:
//
// $Id: DYNAMOFile.C,v 1.7.28.1 2007/03/25 22:00:20 oliver Exp $
//

#include "DYNAMOFile.h"
#include <BALL/DATATYPE/string.h>
#include <BALL/KERNEL/atomContainer.h>
#include <BALL/KERNEL/system.h>
#include <BALL/KERNEL/molecule.h>
#include <BALL/KERNEL/atom.h>
#include <BALL/KERNEL/bond.h>
#include <BALL/KERNEL/PTE.h>
#include <BALL/COMMON/exception.h>
#include <BALL/KERNEL/forEach.h>

#include "../model/Molecule_Visitor.h"
#include "../model/Atom.h"
#include "../model/Molecule.h"
#include "../model/Atom_Element.h"
#include "../model/ModelException.h"
#include "MoleculeReader.h"
#include "BALL_PDB_Types_Adapter.h"

using namespace std;


#include <time.h>
#include <stdlib.h>


/**
 * Returns a speed between -0.5 and 0.5
 */
float get_rand_v()
{
    static bool init = false;
    if (!init)
    {
        srand(time(NULL));
        init = true;
    }

    return ((rand() % 100) / 100.0f) - 0.5;
}

/*  WALLS snippet
"      <!-- DIMENSION X -->\n"
"      <Local Type='Wall' Name='x-min' Elasticity='1' Render='0' Range='All'>\n"
"        <Norm x='1' y='0' z='0'/>\n"
"        <Origin x='-10' y='0' z='0'/>\n"
"      </Local>\n"
"      <Local Type='Wall' Name='x-max' Elasticity='1' Render='0' Range='All'>\n"
"        <Norm x='-1' y='0' z='0'/>\n"
"        <Origin x='10' y='0' z='0'/>\n"
"      </Local>\n"
"      <!-- DIMENSION Y -->\n"
"      <Local Type='Wall' Name='y-min' Elasticity='1' Render='0' Range='All'>\n"
"        <Norm x='0' y='1' z='0'/>\n"
"        <Origin x='0' y='-10' z='0'/>\n"
"      </Local>\n"
"     <Local Type='Wall' Name='y-max' Elasticity='1' Render='0' Range='All'>\n"
"       <Norm x='0' y='-1' z='0'/>\n"
"        <Origin x='0' y='10' z='0'/>\n"
"      </Local>\n"
"      <!-- DIMENSION Z -->\n"
"      <Local Type='Wall' Name='z-min' Elasticity='1' Render='0' Range='All'>\n"
"        <Norm x='0' y='0' z='1'/>\n"
"        <Origin x='00' y='0' z='-10'/>\n"
"      </Local>\n"
"      <Local Type='Wall' Name='z-max' Elasticity='1' Render='0' Range='All'>\n"
"        <Norm x='0' y='0' z='-1'/>\n"
"        <Origin x='0' y='0' z='10'/>\n"
"      </Local>\n"
*/

namespace BALL
{

DYNAMOFile::DYNAMOFile()
    : GenericMolFile()
{
}

DYNAMOFile::DYNAMOFile(const String& name, File::OpenMode open_mode)
throw(Exception::FileNotFound)
    : GenericMolFile()
{
    GenericMolFile::open(name, open_mode);
}

DYNAMOFile::DYNAMOFile(const DYNAMOFile& file)
throw(Exception::FileNotFound)
    : GenericMolFile(file)
{
}

DYNAMOFile::~DYNAMOFile()

{
}

bool DYNAMOFile::write(const AtomContainer& ac)
throw(File::CannotWrite)
{
    if (!isOpen() || getOpenMode() != std::ios::out)
    {
        throw(File::CannotWrite(__FILE__, __LINE__, name_));
    }

    struct : Molecule_Visitor
    {
        std::stringstream hard_sphere_interactions;
        std::stringstream covalent_interactions;
        std::stringstream atom_definitions;

        void on_type_found(const ::Atom_Element* type)
        {
            hard_sphere_interactions <<
                                     "<Interaction Type='HardSphere'"
                                     " Diameter='" << type->radius << "' Elasticity='1'"
                                     " Name='Elem" << type->symbol << "' Range='2All'/>" << endl;
        }

        void on_atom_found(const ::Atom* atom, const ::Atom_Element* type)
        {
            covalent_interactions <<
                                  "<Species"
                                  " Mass='" << type->weight << "'"
                                  " Name='Atom" << atom->id << "'"
                                  " IntName='Elem" << type->symbol << "'"
                                  " Type='Point' Range='Ranged'"
                                  " Start='" << atom->id << "' End='" << atom->id << "'/>" << endl;

            atom_definitions <<
                             "<Pt ID='" << atom->id << "'>"
                             "<P x='" << atom->x << "'"
                             "   y='" << atom->y << "'"
                             "   z='" << atom->z << "' />"
                             "<V x='" << get_rand_v() << "' "
                             "   y='" << get_rand_v() << "' "
                             "   z='" << get_rand_v() << "'/>"
                             "</Pt>" << endl;
        }
    } visitor ;

    BALL_PDB_Types_Adapter pdb(ac);
    MoleculeReader mr;
    ::Molecule molecule;
    // Read the PDB in memory
		try {
    	mr.read(pdb, &molecule);
		}catch(const ModelException &ex){
			std::cout << ex.what() << "\n";
			throw;
		}
    // Translate what was just read to DYNAMO things
    molecule.iterate_atoms_by_type(visitor);

    AtomConstIterator atom_it;
    AtomBondConstIterator bond_it;

    // Compute a "bulk" bond interaction
    double bulk_dmin = +1E+3, bulk_dmax = -1E+3;
    BALL_FOREACH_BOND(ac, atom_it, bond_it)
    {
        // TODO:
        // Consider to use:
        // mili::minimize(bulk_dmin, distance);
        // mili::maximize(bulk_dmax, distance);
        const double distance = bond_it->getFirstAtom()->getDistance(*bond_it->getSecondAtom());
        if (distance < bulk_dmin)
            bulk_dmin = distance;
        else if (distance > bulk_dmax)
            bulk_dmax = distance;
    }
    bulk_dmin *= 0.95;
    bulk_dmax *= 1.05;
    const double bulk_lambda = bulk_dmax / bulk_dmin;

    std::stringstream text;

    text << "<?xml version='1.0'?>\n"
         "<DYNAMOconfig version='1.2.0'>\n"
         "  <Simulation>\n"
         "    <Trajectory Coll='0' nCollPrint='50000'/>\n"
         "    <Ensemble Type='NVE'/>\n"
         "    <Scheduler Type='Dumb'>\n"
         "      <Sorter Type='BoundedPQ'/>\n"
         "    </Scheduler>\n"
         "    <History>Generated by FDMDE's pdb2dynamo.</History>\n"
         "  </Simulation>\n"
         "  <Dynamics>\n"
         "    <Aspect_Ratio x='1' y='1' z='1'/>\n"
         "    <Units Type='SW' BoxLength='1.00000000000000e+02' Energy='1.00000000000000e+00'/>\n"
         "    <BC Boundary='None'/>\n"
         "    <Genus>\n"
         <<
         // "      <Species Mass='1.00000000000000e+00' Name='Bulk' IntName='Bulk' Type='Point' Range='All'/>\n"
         visitor.covalent_interactions.str()
         <<
         "    </Genus>\n"
         "    <Topology/>\n"
         "    <SystemEvents/>\n"
         "    <Globals>\n"
         "      <Global Type='Cells' Lambda='9.00000000000000e-01' Name='SchedulerNBList'/>\n"
         "    </Globals>\n"
         "    <Locals>\n"
         "    </Locals>\n"
         "    <Interactions>\n"
         ;


    // bulk bonds
    text << "      <Interaction Type='SquareBond' Diameter='" << bulk_dmin << "' Lambda='" << bulk_lambda << "' Elasticity='1' Name='Bulk' Range='List'>\n";

    // write all bonds

    BALL_FOREACH_BOND(ac, atom_it, bond_it)
    {
        text << "        <RangePair ID1='"
             << bond_it->getFirstAtom()->getIndex() << "' ID2='"
             << bond_it->getSecondAtom()->getIndex() << "'/>\n";
    }


    text <<
         "      </Interaction>"
         "      <Interaction Type='HardSphere' Diameter='0.5' Elasticity='1.00000000000000e+00' Name='Bulk' Range='2All'/>\n"
         << visitor.hard_sphere_interactions.str() <<
         "    </Interactions>\n"
         "    <Liouvillean Type='Newtonian'/>\n"
         "  </Dynamics>\n"
         "  <ParticleData N='" << (Size)ac.countAtoms() << "' AttachedBinary='N'>\n";


    // write the atom id and the coordinates
    text << visitor.atom_definitions.str();
    text << "  </ParticleData>\n" "</DYNAMOconfig>\n";

    // write all
    (File&)(*this) << text.str();

    return true;
}
bool DYNAMOFile::write(const System& system)
throw(File::CannotWrite)
{
    return write((const AtomContainer&)system);
}

bool DYNAMOFile::write(const Molecule& mol)
throw(File::CannotWrite)
{
    return write((const AtomContainer&)mol);
}

bool DYNAMOFile::read(System& system)
throw(Exception::ParseError)
{
    // remove old rubbish from the system
    system.destroy();

    Molecule* molecule = read();

    if (!molecule)
        return false;

    system.insert(*molecule);
    system.setName(molecule->getName());

    return true;
}

Molecule* DYNAMOFile::read()
throw(Exception::ParseError)
{
    // remember the line number for error messages
    Size number_of_lines = 0;
    String line;

    // prepare a buffer
    const Size BUF_SIZE = 1024;
    char buffer[BUF_SIZE];
    Size number_of_atoms = 0;

    comment_.clear();

    // read the number of atoms (first line)
    if (getline(buffer, BUF_SIZE))
    {
        line.set(buffer);
        line.trim();
        number_of_lines++;
        if (line.countFields() >= 1)
        {
            // retrieve the number of atoms
            number_of_atoms = line.getField(0).toUnsignedInt();
            if (line.countFields() > 1)
            {
                comment_ = line.getField(1);
            }
        }
        else
        {
            // if the first line cannot be read correctly, abort immediately
            Log.error() << "DYNAMOFile::read: illegal header line in XYZ file " << getName() << endl;
            return NULL;
        }
    }
    else
    {
        // we could not read the first line. Abort
        // if the first line cannot be read correctly, abort immediately
        Log.error() << "DYNAMOFile::read: illegal header line in XYZ file " << getName() << endl;
        return NULL;
    }

    // ...create a molecule to hold the atoms, and start reading...
    Molecule* mol = new Molecule;

    // second line: comment -> name of the system
    if (comment_ == "")
    {
        getline(buffer, BUF_SIZE);
        mol->setName(buffer);
        comment_ = buffer;
    }

    Position start = 0;

    HashMap<Position, Atom*> pos_to_atom;

    bool ok = true;

    bool modern_type = 0;

    try
    {
        while (getline(buffer, BUF_SIZE) && (number_of_lines < (number_of_atoms + 2)))
        {
            // read an atom
            number_of_lines++;
            line.set(buffer);

            vector<String> fields;
            Size nr_fields = line.split(fields);
            if (nr_fields < 4)
            {
                Log.error() << "DYNAMOFile: Not enought fields per line!" << std::endl;
                ok = false;
                break;
            }

            if (number_of_lines < 3)
            {
                if (nr_fields > 4)
                {
                    start = 1;
                    modern_type = true;
                }
            }

            // create the atom, insert it into the molecule
            Atom* atom = new Atom;
            mol->insert(*atom);

            // determine the element
            String elementname = line.getField(start);
            if (modern_type)
            {
                atom->setName(elementname);
                elementname = elementname[0];
            }

            Element& element = PTE[elementname];
            if (element == Element::UNKNOWN)
            {
                Log.error() << "DYNAMOFile::read: unknown element " << elementname
                            << " in line " << number_of_lines << " of " << getName()
                            << endl;
            }

            // assign the element and the atom position
            atom->setElement(element);
            atom->setPosition(Vector3(fields[start + 1].toFloat(),
                                      fields[start + 2].toFloat(),
                                      fields[start + 3].toFloat()));

            if (modern_type)
            {
                Position nr = fields[0].toUnsignedInt();
                pos_to_atom[nr] = atom;

                for (Position p = 5; p < nr_fields; p ++)
                {
                    Position partner = fields[p].toUnsignedInt();
                    if (partner < nr)
                    {
                        if (!pos_to_atom.has(partner))
                        {
                            Log.error() << "Could not create bond!" << std::endl;
                            ok = false;
                            break;
                        }

                        atom->createBond(*pos_to_atom[partner]);
                    }
                }
            }
        }
    }
    catch (Exception::GeneralException& e)
    {
        ok = false;
        Log.error() <<  e << std::endl;
    }
    catch (...)
    {
        ok = false;
    }

    if (number_of_lines != number_of_atoms + 1)
    {
        ok = false;
        Log.error() << "Could not read all atoms" << std::endl;
    }

    if (!ok)
    {
        Log.error() << "DYNAMOFile: Aborting, could not parse line!" << std::endl;
        delete mol;
        return NULL;
    }

    return mol;
}

} // namespace BALL
